#!/usr/bin/perl
# Sep 29, 2009 Overlap analysis

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
use Math::Complex;

my $search_width = 300;


if(@ARGV != 2 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [significant region bed file]\n";
}

my %chr_conv = (
	'I'	=>	1,
	'II'	=>	2,
	'III'	=>	3,
	'chr1' => 1,
	'chr2' => 2,
	'chr3' => 3,
);



my %opt;
getopts("i:", \%opt);
my $FILE_target = $opt{i};
#chr1    1478870 1479080 1479080 1259.52 4.296   6.678044e-73    7.26641792070485e-73
#chr1    1479620 1480790 1479860 1380.84 4.553   0       0
#chr1    1484570 1485470 1485020 14211.17        103.211 0       0
#chr1    1526690 1527590 1527140 23857.53        192.027 0       0

my $FILE_tRNA = '/home/hideki/Data/S.Pombe_seq/pombase_ASM294v1.18/contig/tRNA.bed';
my $FILE_5S = '/home/hideki/Data/S.Pombe_seq/pombase_ASM294v1.18/contig/5SrRNA.bed';
my $FILE_LTR = '/home/hideki/Data/S.Pombe_seq/pombase_ASM294v1.18/contig/LTR.bed';
my $FILE_ribosomal = '/home/hideki/Data/S.Pombe_seq/pombase_ASM294v1.18/contig/ribosomal_gene.bed';
#chr1    5064305 5066231 wtf1    0       +
#chr2    593188  593575  wtf2    0       -
#chr3    219185  220283  wtf3    0       -
#chr3    221199  222701  wtf4    0       -
#chr3    243855  244969  wtf5    0       +

my %all;
my %total_number;
my %R_target = &readBed($FILE_target, "target");
my %R_tRNA= &readBed($FILE_tRNA, "tRNA");
my %R_5S = &readBed($FILE_5S, "5SrRNA");
my %R_LTR = &readBed($FILE_LTR, "LTR");
my %R_ribosomal = &readBed($FILE_ribosomal, "ribosomal_gene");
sub readBed{
	my ($file, $title) = @_;
	my %hash;
	my $fh = IO::File->new($file) or die "cannot open $file: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/^track/){
			next;
		}
		my ($chr, $start, $end, @other) = split /\t/;

		unless(exists $chr_conv{$chr}){
			next;
		}

		$total_number{$title}++;

		$all{$chr}{$start - 150}++;
		$all{$chr}{$end + 149}++;
		$hash{$chr}{$start - 150} = 0;	#start
		$hash{$chr}{$end + 149} = 1;		#end
		$hash{'middle'}{$chr}{$start - 150} = int(($start + $end) / 2);	# middleを登録
	}
	$fh->close();
	return %hash;
}


my %peaks_target;
my $summary_tRNA = &analysis(\%R_tRNA, "tRNA");
my $summary_5S = &analysis(\%R_5S, "5SrRNA");
my $summary_LTR = &analysis(\%R_LTR, "LTR");
my $summary_ribosomal = &analysis(\%R_ribosomal, "ribosomal_gene");

#---------------------------------------
# merge peaks
#---------------------------------------
sub analysis{
	my ($Rref, $title) = @_;
	my $flag_inA = 0;
	my $flag_inB = 0;
	my $flag_common = 0;

	my $S_target = 0;
	my $S_basic = 0;
	my $NUM_basic = 0;
	my $TOTAL_basic = 0;

	foreach my $chr(sort {$chr_conv{$a} <=> $chr_conv{$b} } keys %all){
		foreach my $key (sort {$a <=> $b} keys %{$all{$chr}}){
			if(exists $R_target{$chr}{$key} and $R_target{$chr}{$key} == 0){
				$flag_inA = 1;
				$S_target = $key;
				if($flag_inB == 1){
					$flag_common = 1;
				}
			}
			if(exists $Rref->{$chr}{$key} and $Rref->{$chr}{$key} == 0){
				$flag_inB = 1;
				$S_basic = $key;
				if($flag_inA == 1){
					$flag_common = 1;
				}
				$TOTAL_basic++;
			}
			if(exists $R_target{$chr}{$key} and $R_target{$chr}{$key} == 1){
				$flag_inA = 0;
				if($flag_common == 1){
					if($flag_inB == 0){
						# もう１方はすでに終わっている場合はcommon flagを消す
						$flag_common = 0;
					}
					my $distance = abs($Rref->{'middle'}{$chr}{$S_basic} - $R_target{'middle'}{$chr}{$S_target});;
					push @{$peaks_target{$chr}{$S_target}}, "$title\t$distance";
				}
			}
			if(exists $Rref->{$chr}{$key} and $Rref->{$chr}{$key} == 1){
				$flag_inB = 0;
				if($flag_common == 1){
					if($flag_inA == 0){
						# もう１方はすでに終わっている場合はcommon flagを消す
						$flag_common = 0;
					}
					$NUM_basic++;
				}
			}
		}
	}
	my $summary = sprintf "%d / %d = %.2f%%", $NUM_basic, $TOTAL_basic, ($NUM_basic / $TOTAL_basic) * 100;
	return $summary;
}


#---------------------------------------
# １つのピークが２つ以上に帰属している場合には最近接のものを探す(修正：重複の数を数える)
#---------------------------------------
#my %numbers;
my %allpeaks;
my %shared_peaks;
foreach my $chr(keys %peaks_target){
	foreach my $start(keys %{$peaks_target{$chr}}){
		my $MIN = 9999999;
		my $MIN_target = '';

		my @LIST = @{$peaks_target{$chr}{$start}};
		if(@LIST == 1){
			($MIN_target, $MIN) = split /\t/, $LIST[0];
			$allpeaks{$MIN_target}++;
		}else{
			foreach my $key(@LIST){
				my ($target, $dis) = split /\t/, $key;
				$allpeaks{$target}++;
				$shared_peaks{$target}++;
				if($dis < $MIN){
					$MIN = $dis;
					$MIN_target = $target;
				}
			}
		}

#		$numbers{$MIN_target}++;
	}
}

#---------------------------------------
# output result
#---------------------------------------
my $Total_sig_peaks = $total_number{"target"};
print "Peak summary\n";
my $nLTR = exists $allpeaks{'LTR'} ? $allpeaks{'LTR'} : 0;
my $nTRNA = exists $allpeaks{'tRNA'} ? $allpeaks{'tRNA'} :0;
my $n5SrRNA = exists $allpeaks{'5SrRNA'} ? $allpeaks{'5SrRNA'} : 0;
my $nRibosomal = exists $allpeaks{'ribosomal_gene'} ? $allpeaks{'ribosomal_gene'} : 0;

my $sLTR = exists $shared_peaks{'LTR'} ? $shared_peaks{'LTR'} : 0;
my $sTRNA = exists $shared_peaks{'tRNA'} ? $shared_peaks{'tRNA'} :0;
my $s5SrRNA = exists $shared_peaks{'5SrRNA'} ? $shared_peaks{'5SrRNA'} : 0;
my $sRibosomal = exists $shared_peaks{'ribosomal_gene'} ? $shared_peaks{'ribosomal_gene'} : 0;

printf "LTR\t%d (common: %d)\n", $nLTR, $sLTR;
printf "tRNA\t%d (common: %d)\n", $nTRNA, $sTRNA;
printf "5SrRNA\t%d (common: %d)\n", $n5SrRNA, $s5SrRNA;
printf "ribosomal gene\t%d (common: %d)\n", $nRibosomal, $sRibosomal;
printf "others\t%d\n", $Total_sig_peaks - $nLTR - $nTRNA - $n5SrRNA - $nRibosomal;
printf "Total\t%d\n", $Total_sig_peaks;

print "\nCoverage\n";
printf "LTR\t%s\n", $summary_LTR;
printf "tRNA\t%s\n", $summary_tRNA;
printf "5SrRNA\t%s\n", $summary_5S;
printf "ribosomal gene\t%s\n", $summary_ribosomal;






